namespace CSharpKit.Numerics;

using System;

/// <summary>
/// 
/// </summary>
public class Matrix2D
{
    /// <summary>
    /// 构造函数
    /// </summary>
    /// <param name="rows"></param>
    /// <param name="columns"></param>
    /// <exception cref="ArgumentException"></exception>
    public Matrix2D(int rows, int columns)
    {
        if (rows <= 0 || columns <= 0)
            throw new ArgumentException("Dimensions must be positive");

        Rows = rows;
        Columns = columns;
        _data = new double[rows, columns];
    }

    private double[,] _data;
    /// <summary>
    /// 
    /// </summary> <summary>
    /// 
    /// </summary>
    /// <value></value>
    public int Rows { get; }
    /// <summary>
    /// 
    /// </summary> <summary>
    /// 
    /// </summary>
    /// <value></value>
    public int Columns { get; }

    // 索引器
    /// <summary>
    /// 
    /// </summary>
    /// <value></value>
    public double this[int row, int col]
    {
        get => _data[row, col];
        set => _data[row, col] = value;
    }

    // 单位矩阵
    /// <summary>
    /// 
    /// </summary>
    /// <param name="size"></param>
    /// <returns></returns>
    public static Matrix2D Identity(int size)
    {
        var result = new Matrix2D(size, size);
        for (int i = 0; i < size; i++)
        {
            result[i, i] = 1.0;
        }
        return result;
    }

    // 零矩阵
    /// <summary>
    /// 
    /// </summary>
    /// <param name="rows"></param>
    /// <param name="columns"></param>
    /// <returns></returns>
    public static Matrix2D Zero(int rows, int columns)
    {
        return new Matrix2D(rows, columns);
    }

    // 随机矩阵
    /// <summary>
    /// 
    /// </summary>
    /// <param name="rows"></param>
    /// <param name="columns"></param>
    /// <param name="minValue"></param>
    /// <param name="maxValue"></param>
    /// <returns></returns>
    public static Matrix2D Random(int rows, int columns, double minValue = 0.0, double maxValue = 1.0)
    {
        var result = new Matrix2D(rows, columns);
        var random = new Random();
        for (int i = 0; i < rows; i++)
        {
            for (int j = 0; j < columns; j++)
            {
                result[i, j] = minValue + (maxValue - minValue) * random.NextDouble();
            }
        }
        return result;
    }

    // 矩阵转置
    /// <summary>
    /// 
    /// </summary>
    /// <returns></returns>
    public Matrix2D Transpose()
    {
        var result = new Matrix2D(Columns, Rows);
        for (int i = 0; i < Rows; i++)
        {
            for (int j = 0; j < Columns; j++)
            {
                result[j, i] = this[i, j];
            }
        }
        return result;
    }

    // 矩阵加法
    /// <summary>
    /// 
    /// </summary>
    /// <param name="a"></param>
    /// <param name="b"></param>
    /// <returns></returns>
    public static Matrix2D operator +(Matrix2D a, Matrix2D b)
    {
        if (a.Rows != b.Rows || a.Columns != b.Columns)
            throw new ArgumentException("Matrix dimensions must match for addition");

        var result = new Matrix2D(a.Rows, a.Columns);
        for (int i = 0; i < a.Rows; i++)
        {
            for (int j = 0; j < a.Columns; j++)
            {
                result[i, j] = a[i, j] + b[i, j];
            }
        }
        return result;
    }

    // 矩阵减法
    /// <summary>
    /// 
    /// </summary>
    /// <param name="a"></param>
    /// <param name="b"></param>
    /// <returns></returns>
    /// <exception cref="ArgumentException"></exception>
    public static Matrix2D operator -(Matrix2D a, Matrix2D b)
    {
        if (a.Rows != b.Rows || a.Columns != b.Columns)
            throw new ArgumentException("Matrix dimensions must match for subtraction");

        var result = new Matrix2D(a.Rows, a.Columns);
        for (int i = 0; i < a.Rows; i++)
        {
            for (int j = 0; j < a.Columns; j++)
            {
                result[i, j] = a[i, j] - b[i, j];
            }
        }
        return result;
    }

    // 矩阵乘法
    /// <summary>
    /// 
    /// </summary>
    /// <param name="a"></param>
    /// <param name="b"></param>
    /// <returns></returns>
    /// <exception cref="ArgumentException"></exception>
    public static Matrix2D operator *(Matrix2D a, Matrix2D b)
    {
        if (a.Columns != b.Rows)
            throw new ArgumentException("Number of columns in first matrix must match number of rows in second matrix");

        var result = new Matrix2D(a.Rows, b.Columns);
        for (int i = 0; i < a.Rows; i++)
        {
            for (int j = 0; j < b.Columns; j++)
            {
                double sum = 0.0;
                for (int k = 0; k < a.Columns; k++)
                {
                    sum += a[i, k] * b[k, j];
                }
                result[i, j] = sum;
            }
        }
        return result;
    }

    // 标量乘法
    /// <summary>
    /// 
    /// </summary>
    /// <param name="scalar"></param>
    /// <param name="matrix"></param>
    /// <returns></returns>
    public static Matrix2D operator *(double scalar, Matrix2D matrix)
    {
        var result = new Matrix2D(matrix.Rows, matrix.Columns);
        for (int i = 0; i < matrix.Rows; i++)
        {
            for (int j = 0; j < matrix.Columns; j++)
            {
                result[i, j] = scalar * matrix[i, j];
            }
        }
        return result;
    }

    // 标量乘法（交换律）
    /// <summary>
    /// 
    /// </summary>
    /// <param name="matrix"></param>
    /// <param name="scalar"></param>
    /// <returns></returns>
    public static Matrix2D operator *(Matrix2D matrix, double scalar)
    {
        return scalar * matrix;
    }

    // 矩阵求逆（使用高斯-约当消元法）
    /// <summary>
    /// 
    /// </summary>
    /// <returns></returns>
    /// <exception cref="InvalidOperationException"></exception> <summary>
    /// 
    /// </summary>
    /// <returns></returns>
    public Matrix2D Inverse()
    {
        if (Rows != Columns)
            throw new InvalidOperationException("Only square matrices can be inverted");

        int n = Rows;
        var result = new Matrix2D(n, n);
        var temp = new double[n, 2 * n];

        // 创建增广矩阵 [A|I]
        for (int i = 0; i < n; i++)
        {
            for (int j = 0; j < n; j++)
            {
                temp[i, j] = this[i, j];
                temp[i, j + n] = (i == j) ? 1.0 : 0.0;
            }
        }

        // 高斯消元
        for (int k = 0; k < n; k++)
        {
            // 部分主元选择
            double max = Math.Abs(temp[k, k]);
            int row = k;
            for (int i = k + 1; i < n; i++)
            {
                if (Math.Abs(temp[i, k]) > max)
                {
                    max = Math.Abs(temp[i, k]);
                    row = i;
                }
            }

            // 交换行
            if (row != k)
            {
                for (int j = 0; j < 2 * n; j++)
                {
                    (temp[k, j], temp[row, j]) = (temp[row, j], temp[k, j]);
                }
            }

            // 奇异矩阵检查
            if (Math.Abs(temp[k, k]) < 1e-12)
                throw new InvalidOperationException("Matrix is singular and cannot be inverted");

            // 归一化当前行
            double pivot = temp[k, k];
            for (int j = 0; j < 2 * n; j++)
            {
                temp[k, j] /= pivot;
            }

            // 消去其他行
            for (int i = 0; i < n; i++)
            {
                if (i != k)
                {
                    double factor = temp[i, k];
                    for (int j = 0; j < 2 * n; j++)
                    {
                        temp[i, j] -= factor * temp[k, j];
                    }
                }
            }
        }

        // 提取逆矩阵
        for (int i = 0; i < n; i++)
        {
            for (int j = 0; j < n; j++)
            {
                result[i, j] = temp[i, j + n];
            }
        }

        return result;
    }

    // 行列式计算（使用递归方法）
    /// <summary>
    /// 
    /// </summary>
    /// <returns></returns>
    public double Determinant()
    {
        if (Rows != Columns)
            throw new InvalidOperationException("Only square matrices have determinants");

        int n = Rows;

        // 1x1矩阵
        if (n == 1) return this[0, 0];

        // 2x2矩阵
        if (n == 2) return this[0, 0] * this[1, 1] - this[0, 1] * this[1, 0];

        double det = 0.0;
        for (int j = 0; j < n; j++)
        {
            det += (j % 2 == 0 ? 1 : -1) * this[0, j] * Minor(0, j).Determinant();
        }
        return det;
    }

    // 获取子矩阵（用于计算余子式）
    /// <summary>
    /// 
    /// </summary>
    /// <param name="rowToRemove"></param>
    /// <param name="colToRemove"></param>
    /// <returns></returns>
    private Matrix2D Minor(int rowToRemove, int colToRemove)
    {
        var result = new Matrix2D(Rows - 1, Columns - 1);
        int r = 0;
        for (int i = 0; i < Rows; i++)
        {
            if (i == rowToRemove) continue;

            int c = 0;
            for (int j = 0; j < Columns; j++)
            {
                if (j == colToRemove) continue;

                result[r, c] = this[i, j];
                c++;
            }
            r++;
        }
        return result;
    }

    // 矩阵转换为字符串
    /// <summary>
    /// 
    /// </summary>
    /// <returns></returns>
    public override string ToString()
    {
        var sb = new System.Text.StringBuilder();
        for (int i = 0; i < Rows; i++)
        {
            sb.Append("[ ");
            for (int j = 0; j < Columns; j++)
            {
                sb.Append(this[i, j].ToString("F4")).Append(" ");
            }
            sb.AppendLine("]");
        }
        return sb.ToString();
    }

    // {{END}}
}